home *** CD-ROM | disk | FTP | other *** search
/ CD-ROM Today - The Disc! 5 / CD-ROM Today - The Disc (Issue 5)(November 1994).ISO / mac / Mac shareware / Education / RLaB / examples / house.r < prev    next >
Text File  |  1994-09-21  |  2KB  |  118 lines
  1. //
  2. // From MATRIX Computations, G.H. Golub, C.F. Van Loan
  3. //
  4.  
  5. //
  6. // house_v(): Given an N-vector V, generate an n-vector V
  7. // with V[1] = 1, such that (eye(n,n) - 2*(V*V')/(V'*V))*X
  8. // is zero in all but the 1st component.
  9. //
  10.  
  11. static (sign)
  12.  
  13. house_v = function(v)
  14. {
  15.   local(v, b, n, u)
  16.   
  17.   n = max( size(v) );
  18.   u = norm(v, "2");
  19.   if(u != 0)
  20.   {
  21.     b = v[1] + sign(v[1])*u;
  22.     if(n > 1) 
  23.     {
  24.       v[2:n] = v[2:n]/b;
  25.     }
  26.   }
  27.   v[1] = 1;
  28.   return v;
  29. };
  30.       
  31. //
  32. // house_row(): Given an MxN matrix A and a non-zero M-vector V
  33. // with V[1] = 1, the following algorithm overwrites A with 
  34. // P*A, where P = eye(m,m) - 2*(V*V')/(V'*V)
  35. //
  36.  
  37. house_row = function(A, v)
  38. {
  39.   local(A, b, w)
  40.   
  41.   b = -2/(v'*v);
  42.   w = b*A'*v;
  43.   A = A + v*w';
  44.   return A;
  45. };
  46.  
  47. //
  48. // house_col(): Given an MxN matrix A, and an N-vector V, 
  49. // with V[1] = 1, the following algorithm overwrites A with A*P
  50. // where P = eye(N,N) - 2*(V*V')/(V'*V)
  51. //
  52.  
  53. house_col = function(A, v)
  54. {
  55.   local(A, b, w)
  56.   
  57.   b = -2/(v'*v);
  58.   w = b*A*v;
  59.   a = A + w*v';
  60.   
  61.   return A;
  62. };
  63.  
  64. //
  65. // Given A, with M >= N, the following function finds Householder
  66. // matrices H1,...Hn, such that if Q = H1*...Hn, then Q'*A = R is
  67. // upper triangular.
  68.  
  69. // House_qr returns a MxN matrix, with the upper triangular part 
  70. // containing [R]
  71.  
  72. house_qr = function ( A )
  73. {
  74.   local (A, j, n, m, q, v)
  75.   
  76.   m = A.nr; n = A.nc;
  77.   v = zeros(m,1);
  78.   q = eye (m, m);
  79.  
  80.   for(j in 1:n)
  81.   {
  82.     // Generate the Householder vector
  83.     v[j:m] = house_v( A[j:m;j] );
  84.  
  85.     // Apply the Householder reflector to A
  86.     A[j:m;j:n] = house_row( A[j:m;j:n], v[j:m] );
  87.  
  88.     // Create Q
  89.     if(j < m) 
  90.     {
  91.       q = P ([ zeros (j-1,1); 1; v[j+1:m] ]) * q;
  92.     }
  93.   }
  94.   return << q = q'; r = A >>;
  95. };
  96.  
  97. //
  98. //  Generate P matrix
  99. //
  100.  
  101. P = function ( V )
  102. {
  103.   local( m )
  104.   
  105.   m = max( size(V) );
  106.   return [ eye(m,m) - 2*(V*V')./(V'*V) ];
  107. };
  108.  
  109. sign = function ( X )
  110. {
  111.   if (X >= 0)
  112.   { 
  113.     return 1; 
  114.   else 
  115.     return -1;
  116.   }
  117. };
  118.